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By detailed Molecular Dynamics and Monte Carlo simulations we show that granular materials 
at rest can be described as thermodynamics systems. First we show that granular packs can be 
characterized by few parameters, as much as fluids or solids. Then, in a second independent step, 
we demonstrate that these states can be described in terms of equilibrium distributions which 
coincide with the Statistical Mechanics of powders first proposed by Edwards. We also derive 
the system equation of state as a function of the "configurational temperature" , its new intensive 
thermodynamic parameter. 



Granular materials, such as sand or powders, for many 
respect are similar to fluids or solids [l|, even though in 
absence of external drive they rapidly come rest, due to 
strong dissipation and negligible thermal energy scales 
(they are non-thermal systems ), in disordered states very 
similar to glasses [MSSIMIEli^- As standard thermo- 
dynamics is not applicable to describe them, it is natu- 
ral to ask whether we can even refer to granular packs 
as "states". The problem of finding the correct theo- 
retical framework where to describe granular media is in 
fact of deeprelevance to civil engineering, geophysics and 
physics llaili 

Edwards proposed a thermodynamic description 
for static granular media, which was partially investi- 
gated by recent experiments 0, Q ■ These exper- 
iments have established that a granular system subject 
to a tapping dynamics, such as subsequent mechanical 
oscillations of the container, may loose memory of its ini- 
tial state and reach a stationary state of volume fraction 
only dependent on the tapping intensity, a precondition 
for a statistical mechanics description of static granular 
material to be possible. The study of out-of equilibrium 
(aging) slowly sheared granular assemblies is also useful 
for the validation of the statistical mechanics of granu- 
lar media 0|, but it is inherently restricted to a small 
range of very high volume fractions, where the system is 
jammed [l^. 

Here we give strong evidences supporting the existence 
of a thermodynamical and statistical mechanical descrip- 
tion of granular media. First we demonstrate, via Molec- 
ular Dynamics (MD) simulations, that granular packs at 
rest are genuine thermodynamic states, as they are char- 
acterized by a small set of parameteres regardless of the 
procedure with which they are generated. Then we show, 
via Monte Carlo (MC) simulations, that these states can 
be described in terms of the equilibrium distribution pro- 
posed by Edwards. The coincidence between time aver- 
ages (MD) and ensemble averages (MC) is a strong evi- 
dence in favour of the statistical mechanics approach to 
granular media. For details of materials and methods, 
see the supplementary materials available online 

MD simulations: time averages - We run Molecular 
Dynamics simulations of = 1600 monodisperse spheri- 
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FIG. 1: (color online) Schematic plot of the investigated 
system. A monodisperse 3D grain pack of spheres immersed 
in a fluid is subject to a dynamics made of flow pulses of 
velocity V and duration tq. Before applying a pulse we wait 
for the pack to settle. The flow pulse dynamics is repeated 
up to reach stationary conditions (see Figl^J. 



cal grains of diameter d = 1cm and mass to = Ig. Grains, 
under gravity, are confined in a box with a square basis 
of length L = 10cm (see Fig.^, with periodic boundary 
conditions in the horizontal directions. The bottom of 
the box is made of other immobile, randomly displaced, 
grains (to prevent crystallization) . Two grains in contact 
interact via a normal and a tangential force. The former 
is given by the spring-dashpot model, while the latter is 
implemented by keeping track of the elastic shear dis- 
placement throughout the lifetime of a contact [iR ll^ . 
The coefficient of restitution is constant, e = 0.8. 

As in a recent experiment the system is immersed 
in a fluid and, starting from a random configuration, it is 
subject to a dynamics made of a sequence of flow pulses 
where the fluid flows through the grains (see Fig.^. In 
a single pulse the flow velocity, directed against gravity, 
is y > for a time tq; then the fluid comes to rest. 
We model the fluid-grain interaction 0, via a vis- 
cous force proportional to the fluid grain relative veloc- 
ity: Ffg = —A{'v — V) where v is the grain and V is 
the fluid velocity. The prefactor A = 7(1 — is 
dependent on the local packing fraction, , in a cube of 
side length 3d around the grain, and the constant d9] is 
7=1 Ns/cm. 

During each pulse, grains are fluidized and then come 
to rest under the effect of gravity. The tapping dynamics, 
therefore, allows for the exploration of the phase space of 
the mechanically stable granular packs. When the system 
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FIG. 2: (color online) Compaction of systems subject to 
a tapping dynamics with to = 0.03s and different values of 
the fluid velocity as indicated, averaged over 32 runs. Stars 
refers to simulations with 4 times more particles, outlining 
the absence of finite size effects. Dashed lines are fits to a 
stretched exponential law, Eq. Q 



is subject to such a tap dynamics, it compactifies until 
it reaches a stationary state where its properties do not 
depend on the dynamics history. Fig. El shows that the 
volume fraction of our system increases by following a 
stretched exponential low. 



= $oo - ($oo - $o) exp {-{t/Tf ) , 



(1) 



in agreement with the experiment by P. Philippe et 
al. |20j. The relaxation time diverges as the tapping 
intensity decreases, indicating the presence of a glassy 
like behavior which will be discussed elsewhere [2l|. As 
the thermodynamics approach to granular media aims to 
describe stationary states, all of measures shown below 
(averaged over 32 runs) are recorded after the application 
of a long sequence of flow pulses, when the system is at 
stationarity. 

We plot in Fig. |3| the stationary values of the volume 
fraction, ^(V,to) (measured in the bulk of the system), 
and its fluctuations, ^(jiiV, Tq), recorded after a sequence 
of such flow pulses of duration tq and velocity V . A(/) is 
by definition the standard deviation of (j) around its av- 
erage value at stationarity. Actually, the volume fraction 
probability distributions is Gaussian 0,^3- decreases 
with V and with tq: the stronger the pulse, i.e., the larger 
1^ or To, the fluffier the pack settled after it. Similarly, 
increases with V and tq. 

Even though, depends on both the parameters of the 
dynamics, V and tq, we show now that such stationary 
states are indeed genuine "thermodynamic states", i.e., 
they can be described, in this system, by one macroscopic 
parameter. Actually, the upper panel of Fig. 0] shows 
that when A0 is parametrically plotted as function of 
(/), the scattered data of Fig. |31 collapse, within numeri- 
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FIG. 3: (color online) The main panel shows the volume 
fraction (j) attained at stationarity (see Fig.j^J measured when 
the pack is settled, as a function of V , for the shown values 
of Tq: the stronger the pulses the fiuffier the pack. The inset 
shows the dependence of the volume fraction fiuctuations on 
V and Tq. 



cal approximation, onto a single master function. This 
is a clear indication that A</) and (f) are in a one-to-one 
correspondence, no matter how the state with packing 
fraction is attained. Our claim is that such a prop- 
erty should be found for any macroscopic observable of 
the system: we checked some of them, including the en- 
ergy and its fluctuations and the coordination number of 
grains. Actually, in Fig. 0] lower panel we show that the 
whole radial distribution function g{r) of a pack is char- 
acterized only by its corresponding value of 0, i.e., states 
attained with different dynamical protocols (y, tq), but 
having the same 0, have the same g{r). From these re- 
sults we derive our first conclusion: at stationarity, we 
can describe the pack with only one parameter, e.g., 0, 
independently of the dynamical protocol. Such a param- 
eter characterizes, thus, the "thermodynamics state" of 
the system. 

MC simulations: ensamble averages - Our second im- 
portant step is to identify the correct Statistical Mechan- 
ics distribution for these states. Under a very strong as- 
sumption (discussed for instance in |0, 0, 

Edwards proposed to use for 
the grains of a powder the standard machinery of Statis- 
tical Mechanics. He suggested, however, to consider a re- 
duced configurational space: the system at rest (i.e., not 
in its "fluidized" regime) is described by a flat ensemble 
average restricted to its blocked configurations (i.e., its 
mechanically stable microstates). Under these hypothe- 
ses 0, 13 Ell the canonical ensemble probability, Pr , to 
find the blocked microstate r, of 6ii6rgy is! 



p. 



oc e 



(2) 



3 





0.30 r 


-e- 0.27 - 


<1 










0.24- 








0.21- 












0.18- 


'o 


> 






0.15- 




A Xq = 0.03 s 
Xq = 0.06 s 
• Tq = 0.09 s 
■ Xq = 0.12 s 
O Ens. Averag. 



^^^^ 



0.57 0.58 0.59 0.6 0.61 

Volume fraction, (|) 




1.6 1.8 2 2.2 

Interparticle distance, r (cm) 

FIG. 4: (color online) Upper panel The A(j) data of Fig. 01 
(same symbols are used) are here plotted as a function of (p: 
they collapse on a single master curve, showing that A<j> is 
in a one-to-one correspondence with <j!>, irrespectively of the 
dynamical protocol used to arrive at </!>. The dotted line is a 
linear best fit. A similarly good collapse is also found (lower 
panel) when we plot the radial distribution function g{r) for 
packs having the same volume fraction (the same packs of 
Fig. 121 are shown with the same symbols). In particular, we 
show g{r) for (f) = 0.615 and (/> = 0.575 (shifted for clarity). 
The empty circles in both panels are the corresponding en- 
semble averages independently calculated from Eq.|51 within 
numerical errors, they scale on the same curves, pointing out 
this Statistical Mechanics measure is an excellent approxima- 
tion of the time averages over the pulse dynamics. 



where the inverse of Pconf is the conjugate parameter 
of the energy, called conf igurational temperature. 

In order to check whether such a Statistical Mechan- 
ics scenario applies, we compared ensemble averages over 
the distribution of Eq.|21with those over the flow tap dy- 
namics. For instance, the average value of over the 
distribution of Eq. |21 is 



{(t))iTco„f) 



X],. (I>r exp{-Er/Tconf) 
J2r'^^P(~^r/Tco7if) 



(3) 



where the sum runs over all blocked microstates, and 
(j)r is the volume fraction of microstate r. We evaluated 
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FIG. 5: (color online) The equation of state of granular ma- 
terials at rest, showing the volume fraction, (j>, as a function 
of the configurational temperature, Tconf- Ensamble averages 
are obtained via the Monte Carlo procedure, while time av- 
erages are obtained from the data of Fig. El via the use of the 
static fluctuation dissipation relation |15| . 



these ensemble averages by use of a Monte Carlo method 
which is an extension of that introduced in Ref. 6] to 
the frictional case (frictional forces are essential to 
assure the stability of granular packs with small vol- 
ume fraction). Fig. 01 shows, as empty circles, the fmic- 
tions (A0)((0)) (resp. {g{r)){{(f>))) in the upper (resp. 
lower) panel. These ensemble averages collapse, to a 
very good approximation, on the same master function 
of the time averaged data from the flow-pulse dynam- 
ics discussed before (notice that there are no adjustable 
parameters). Thus, the present Statistical Mechanics de- 
scription appears to hold, up to the current numerical 
accuracy, at least as a first very good approximation. In- 
terestingly, the off-equilibrium dynamical effective tem- 
perature defined at high volume fractions from dynam- 
ical fluctuation-dissipation relations appears to coin- 
cide with the configurational temperature derived here 
at stationarity 

The function (t>(Tconf), derived from the above ensem- 
ble averages, is the system equation of state. We plot 
4'{Tconf) in Fig. [31 where the data from MD simulations 
are included too by using the data collapse from Fig. ^ 
and i nteg ration of the static fiuctuation-dissipation rela- 
tions lilElElillii. 

Conclusions - Summarizing, in the present MD sim- 
ulations of a non-thermal monodisperse granular system 
under fiow pulses, we find that the stationary configura- 
tions of the system can be fully described by only one 
parameter, e.g., 0, and can be, thus, considered genuine 
"thermodynamic states" . Within our numerical accu- 
racy, we also showed that a Statistical Mechanics based 
on the distribution of Eq.Elis grounded to describe these 
"states". We could derive as well the equation of state, 
4>{Tconf), of the system. 
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These evidences strongly support the existence of a 
fundamental theory of dense granular media and address, 
thus, a variety of important issues in the next future, 
such as response functions in a granular system, mix- 
ing/segregation phenomena, the nature of their jamming 
transition and phase diagram 
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Volume fraction fluctuations in the stationary state 




Volume fraction, O 



Figure A Probability distributions of the volume fraction $ in the steady state reached 
by the system subject to the flow pulses dynamics. The data, obtained by averaging over 
the stationary dynamics in 32 different runs starting from random initial conditions, are 
Gaussian distributed (plain lines). The standard deviations of these distributions are the 
volume fraction fluctuations, A<i>, plotted in Fig. 3. 



Statistical Mechanics approach to granular media 

In the Statistical Mechanics of powders introduced by S. Edwards (S.F. Edwards and R.S.B. 
Oakeshott, Physica A 157, 1080, 1989) it is postulated that the system at rest can be de- 
scribed by suitable ensemble averages over its "mechanically stable" states. Edwards pro- 
posed a method to individuate the probability, Pr, to find the system in its mechanically 
stable state r, under the assumption that these mechanically stable states have the same 
a priori probability to occur. This is the simplest assumption one can imagine: it is the 
assumption of standard Statistical Mechanics (equiprobability of microstates) with the ad- 
ditional constraint of mechanical stability. A possible approach to find Pr is as follow (A. 
Fierro et ai, Europhys. Lett. 59, 642, 2002; Phys. Rev. E 66, 061301, 2002; Europhys. 
Lett. 60 684, 2002). is obtained as the maximum of the entropy, 

S=-J2Prl0gPr 

r 

with the macroscopic constraint, in the case of the canonical ensemble, of fixed system 
energy E = J^PrEr (here Er is the energy of the mechanically stable microstate r). For a 
granular medium at rest in the gravity field, Er is the sum of the gravitational energy, and 
of the interaction energy between grains; 

Er = ^ mgzi + Vij , (4) 

where Zi is the height of grain i (with mass m) and Vij is the interaction (elastic energy) 
between grain i and j in the microstate r: Vij = if grains i and j are not in contact, 
otherwise 

Vij = ^^nl4P + \h\Uij\^, (5) 

where Sij is the overlap between grain i and j, and Uij is their shear displacement. The 
interaction energy Vjj is derived by the normal and the tangential forces acting on the 
contacting grains (we use the model L3 described in L.E. Silbert et al, Phys Rev E 64, 
051302, 2001): /" = fc„4. = ^A-j. 

From Edwards' hypothesis, in analogy to the Gibbs result, you derive that: 

Pr = exp {-PconfEr) 

where f3conf is a Lagrange multiplier, called inverse configurational temperature, enforcing 
the above constraint on the energy: 

Pconf = — qI^'^ ; Sconf = ln f^stable (-£') • 

The configurational temperature is Tconf = Pconf Here instable (-E-) is the number of 
mechanically stable states with energy E. Z is fixed by the normalization condition 
Pr = 1, where the sum is restricted to rigtabie; i-e., to the mechanically stable states. 

We note here that for stiff grains the elastic energy Y^i^j Vij is much smaller than 
the gravitational energy mgzi (the elastic energy is strictly zero in hard sphere systems) . 
Accordingly Er is to a good approximation proportional to the gravitational energy, i.e. to 
the volume of the system, as originally suggested by Edwards. 



Configurational temperature and dynamical tempera- 
ture 

We describe here the two different definitions of temperature for granular systems mentioned 
in our paper. 

Configurational temperature 

As discussed above, the configurational temperature is the inverse of the derivative of 
the configurational entropy with respect to the energy. It is therefore an 'equilibrium' 
temperature, defined for a granular system at stationarity under a given dynamics allowing 
the exploration of mechanically stable states. 

In our main text we have validated Statistical Mechanics approaches to powders by 
showing that, within numerical errors, the granular packs we consider are characterized, 
at stationarity, by a single parameter regardless of the dynamical procedure with which 
they were prepared (they are in a thermodynamic 'state'), and that time averages coincide 
with ensemble averages over the distribution detailed above. This result justifies the use of 
Edwards 'equilibrium' partition function to make analytical calculations, and particularly 
to derive the following equilibrium fluctuation-dissipation relation (Coniglio and Nicodemi 
Physica A 296, 451 (2001)), 



which relates the energy and its fluctuations. This relation can be exploited to experimen- 
tally measure the conflgurational temperature (Nowak et al., Powder. Tech. 94, 79, 1997; 
Knight et al, Phys. Rev. E 57, 1971, 1998; Schroeter and Swinney, Phys. Rev. E 71, 
030301 (R), 2005). 

Dyneunical temperature 

In thermal glassy systems far from stationarity, dynamical off-equilibrium fluctuation- 
dissipation relations hold. Particularly, in the aging dynamics of mean-field glassy models 
in contact with a very small bath temperature, Tbath, generalized out-of equilibrium 
fluctuation-dissipation relations were discovered where the role of the external temperature 
is played by a 'dynamical temperature' Tdyn 7^ ^bath, equal for all slow modes (L.F. 
Cugliandolo and J. Kurchan, Phys. Rev. Lett. 71, 173, 1993). This scenario was later 
extended to aging gramdar materials (for a review sec P. Richard et a/., Nat. Mat. 4, 121, 
2005, and references therein), were the dynamical temperature is defined and measured via 
a dynamical, off-equilibrium fluctuation dissipation relation, 



where r is the position of a grain, / a constant perturbing fleld, and tw the 'waiting' time. 
As the dynamical temperature is defined far from stationarity and, conversely, the configu- 
rational temperature at stationarity, Makse and Kurchan have shown that in granular packs, 
at very high density, the two are equal within numerical errors (Nature 415, 614, 2002). At 
low volume fraction, on the contrary, the dynamical temperature appears to be no longer 
defined, as the granular system is no longer jammed and fiows as soon as an external per- 
turbation is applied, as shown by Potiguar and Makse (European Physical Journal E 19, 




((r(i)-r(i^))')=Tdyn 



S{r{t) - r{t^)) 
Sf 



171, 2006). 



Monte Carlo method to test the Statistical Mechanics 
approach to granular materials at rest 



In order to test the Statistical Mechanics approach to granular media one has to 
compare time averaged data, obtained as explained in the text, with ensemble averaged 
data obtained by sampling Edwards distribution, eq.(2) in the text, 

Pr^ Z-^eM-PconfEr). (6) 

where Er is the energy of the state r, 

= ^m.gzi + ^l^y. (7) 

As the Statistical Mechanics approach to granular media deals with mechanically stable 
states, the phase space of interest (over which eq. ^ must be sampled) is not the usual phase 
space, f2tot{?', "S", w} (here r is the iN vector of grains c.o.m. positions, v their velocities, 
and uj their angular velocities). Instead it is the subset Jlstabic C Q.tot{?^ 0, 0} of all states r 
where the forces and the torques acting on each single grain sum to zero, and grains neither 
translate nor rotate. 

Since the states to be considered are so highly constrained it is difficult to sample the 
distribution of eq. lO via a standard Monte Carlo procedure. For instance, if a state r is 
stable, then there is little chance to transform it in a new mechanically stable state r' via 
the displacement of a single particle. Introducing many particles Monte Carlo moves is also 
useless as the probability of selecting a collective move that transform a mechanically stable 
state into a new mechanically stable state is practically zero. 

A Monte Carlo (MC) method to explore Instable was proposed by H.A. Makse and J. Kurchan, 
Nature 415 614 (2002), which uses the following computational trick. The MC algorithm 
explores the usual phase space, lltot, but in this phase space one introduces an auxiliary 
energy Eaux which measures the degree of 'mechanical instability' of a pack in a microstate 
r. In the present case, we have defined: 

N N 

Eaux = ^'^ad + ^(4;- + 4)1 + i^^i' (8) 

2—1 j^i i—1 

where is the total torque acting on grain i. The first term of the above equation enforce 
the balance of forces on each single grain, and the second the balance of torques. Different 
expressions could be used for Eaux, the important point being that Eaux{r) > Vr, and 
that Eaux (r) = if and only if the state r is mechanically stable. 

Then one samples via a standard Monte Carlo procedure the distribution 

Pr = Z ^ exp{~(3confEr - PauxEaux{r)) , (9) 

where Taux = Paux the so-called 'auxiliary temperature' which controls the equilibrium 
value of the auxiliary energy. By definition of Eaux, in the limit Ta„x ^ we have: 

Pr = hm Pr, hm Eaux = 0. (10) 

Therefore in the limit Taux we sample the distribution of mechanically stable states 
(Eaux = 0) with probability Pr, as desired. This is precisely the limit which is considered 
in by Makse and Kurchan method. 

In our simulations we start by sampling eq. ||^ in the phase space of all granular packs at 
the desired value of Tconf and at Taux > via a Monte Carlo procedure (see below). Then 



we slowly decrease the auxiliary temperature (carefully checking for thermalization) until 
Taux = 0. By repeating this procedure severals times we generate several packs (a total of 
172) at the desired configurational temperature. These packs are then used to compute 
ensemble averages relative to the chosen value of Tconf- 

Implementation of the Monte Carlo method 

In our definition of mechanical stability (cq. (jSJ) we also include tangential forces (neglected 
in H.A. Makse and J. Kurchan, Nature 415 614, 2002). Tangential forces model friction, 
which is essential for our purposes since, under gravity, frictionless stable packs are only 
found for high volume fractions. 

In MD simulations friction has important effects on the dynamics of the system: the 
frictional force between two grains depends on their shear displacement. In the Monte Carlo 
algorithm, however, it is convenient to consider the frictional force between two grains (that 
is their shear displacement) as an independent variable. This leads to the definition of two 
kind of MC moves: standard single particle displacement, and variation of the tangential 
shear displacement. The idea of separating geometrical properties of the pack (particle 
displacement) from the tangential forces (tangential shear variation) has been exploited 
previously in the literature, even though in a different contest and with simpler models (as 
for example in T.Unger, J.Kertesz, and D. E. Wolf, Phys. Rev. Lett. 94, 178001 (2005); S. 
McNamara, R. Garcma-Rojo, and H. Herrmann, Phys. Rev. E 72, 021304 (2005)). Below 
we discuss briefly the two moves. 

Single-particle displacement 

One selects a particle n in position r„ and a random displacement vector A, with |A| < A 
and A dynamically varied in order to obtain an acceptance probability Pacc = 0.5. The 
displacement of particle n from position r„ to position r„ + A induces a variation AEr of 
the energy (eq. I@J), and of the auxihary energy, AEaux (eq. (|H1)), of the system. AEaux is 
due to: 

1. changes of the overlaps 5ni between the displaced particle n and a particle i, and 
therefore of the normal force /^j. In particular contacts may disappear (after the 
displacement \8ni < 0|), or appear (before the displacement \Sni < 0|, after the dis- 
placement \Sni > 0|). 

2. variations of the shear displacements Uni, and therefore of the tangential force f^^. If 
particles i and n are in contact both before and after the displacement of particles 
n, we consider particle n as sliding over particle i, inducing a variation of the shear 
displacement Uni- If the displacement of particle n creates (destroys) a contact, then 
the shear displacement u„i is created (destroyed) accordingly. 

The shear displacement is rescaled by a factor ^|/j"|/A:t|u„i| if the Coulomb condition is 
violated. A move of this kind is accepted with probability exjp{— f3conf AEr — PauxAEaux)- 

Shear-displacements variation 

In this MC move one varies the tangential force between two contacting grains, respecting 
the Coulomb criterion. Two particles n and m are selected. If they are in contact one varies 
their shear displacement Unm of a random amount At, with |At| < a and a dynamically 
varied in order to obtain an acceptance probability pacc = 0.5. At is actually chosen in such 
a way that Unm + A( (and therefore the tangential force) lies in the plane tangent to both 
grains in their point of contact. Before varying by At this latter is rescaled in order to 
satisfy the Coulomb criterion, if necessary. 

The tangential move leads to a variation of both the energy of the system, Er (as the elastic 
energy in tangential interaction between the two particles varies), and of its auxiliary en- 
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ergy, Eaux- As before, the move is accepted with probabihty exp(— /^con/A-Er — Paux^Eaux)- 
Generation of mechanically stable states with the MC algorithm 
0.62 
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Here we show how the volume fraction of the system evolves during the MC procedure. 
First we thermalize the system at the desired value of Tconf and at Taux = 0.04 mgd (upper 
panel). Different values of Taux may be convenient in different runs. Then we slowly decrease 
the auxiliary temperature until the system reaches a mechanically stable state (lower panel) . 
We check that the same results ar e obtained with slower cooling rates. 



